The Kuramoto Model with Time- Varying Parameters 
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We introduce a generalization of the Kuramoto model by explicit consideration of deterministically 
time-varying parameters. The oscillators' natural frequencies and/or couplings are influenced by 
external forces with constant or distributed strengths. We observe a new dynamics of the collective 
rhythms, consisting of the external system superimposed on the autonomous one, a characteristic 
feature of many thermodynamically open systems. This deterministic, stable, continuously time- 
dependent, collective behaviour is fully described, and the external impact is defined in both the 
adiabatic and non-adiabatic limits. 
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Biological examples provided the original motivation 
lying behind the Kuramoto model (KM) of coupled 
phase oscillators [l|. However, neither the original model 
[4I, nor any of its extensions [3||, have incorporated a 
fundamental property of living systems - their inher- 
ent time-variability. Many important characteristics of 
open systems can be missed by not accounting for the 
non-equilibrium dynamics that stems from their time- 
dependent (TD) parameters. Additionally, the applica- 
tion of the KM to many problems would move closer to 
reality by allowing for the natural frequency of each os- 
cillator, or the coupling strengths, to be externally mod- 
ulated by TD forcing, as commonly occurs in living sys- 
tems. Among the numerous collective rhythms traceable 
back to TD parameters, are frequency flows in brain sig- 
nals [4, the modeling of brain dynamics under anaes- 
thesia [5| where anaesthetic strength modulates natural 
frequencies [6j , event-related oscillatory responses of the 
brain [7| , and the dynamics of cardiovascular ageing [8| . 
None these is adequately described by existing models. 

We introduce an external, explicitly time-dependent, 
bounded function x{t), modulating the frequencies or 
couplings of the original model. The influence can also 
come from another [9[ non-constant mean field. In the 
most general case, the strengths of the interactions li are 
distributed according to a probability density function 
h(I), and likewise the distribution g{ijj) of the natural 
frequencies uii. Thus, depending on which parameter is 
influenced two generalized Kuramoto models emerge 

A: e,=io,+hx{t)+K r{t)sm{^-e^), (1) 

B: 9,=uj, + [K + I,x{t)]r{t)sm{^-e,). (2) 

Here, we have introduced a TD complex order parameter 

1 ^ 
zit) = r{t)e^m^-J2e^^^, (3) 

where r(t) and ip{t) ^^^ the TD mean- field amplitude 
and phase respectively. For clarity, their explicit time- 
dependence will henceforth be omitted. 

There has already been much work on coupled oscil- 
lators influenced by noise as a special form of external 



dynamics [10|. Other studies of non-constant collective 
rhythms include asymmetrically-coupled ensembles 111. 



and models with inertia 12|, frequency adaptation 13|, 
alternately-switching 1J| or periodic couplings [l5j, or 
driving by an external periodic force [16|. However, the 
TD mean fields in these models result either from mul- 
tistability or from unstable equilibria. As such, none of 
them can demonstrate the deterministic and stable TD 
dynamics of many real physical, chemical, biological, or 
social systems that can never be completely isolated from 
their surroundings. These systems do not reach equilib- 
rium but, instead, exhibit complex dynamical behavior 
that includes the TD frequencies and couplings. We will 
show that our generalization of the Kuramoto model en- 
compasses this dynamics. 

To analyze the models ([1]), ([2]) we assume the ther- 
modynamic limit A'^ ^ 00, where there is 1:1 corre- 
spondence between the fixed and TD parameters, i.e. 
LJ{uj,I,t) = Lu + Ix{t) for model A, and k{K,I,t) = 
K + Ix{t) for model B. Therefore, the state of the oscil- 
latory system can be described by a continuous prob- 
ability density function p{9,uj,I,t) [17| normalized as 
/p p{9,uj,I,t)d9 = 1. It gives the proportion of oscil- 
lators with phase 9 and at time t, for particular w and 
/. Hence, conservation of the number of oscillators gives 
the continuity equation for every fixed uj and / 

Here, we used the definition (|3]), rewritten using 
W Si sin(6'j — 9i) = Ini{ze~'^'}, so that it now becomes 

^27r /"OO /'C30 

/ p{uj,I,9,t)g{ij)h{I)e''^d9dujdI. (5) 

00 ^ — c« 

Since pico, I,9,t) is real and 27r periodic in 9, it allows a 
Fourier expansion with the Ott and Antonsen ansatz [l8| 
fn{u>,I,t) = [a{uj , I , t)]"' in its coefficients. Thus, 

pi9, u;,I,t)^—{l + {J2 H^, I, i)]"e™' + cc}}. (6) 



n=l 



where c.c. stands for complex conjugate, is a particular 
solution of Eqs. (jU [5]), as long as a{uj, I, t) evolves with 

A: ^ + i[oj + lxit)]a + ^{za^-z*)^0, (7) 



B: 



da 
'dt 



K + Ix{t) ^ 2 



{za^ - z*) = 0, (8) 



for models A and B respectively. The same ansatz im- 
plemented in Eq. ([S]), reduces the order parameter to 



a(a;, /, t)g{uj)h{I)dLodI . 



(9) 



In all further analysis, the natural frequencies follow 
a Lorentizan distribution, and a{ui,I,t) is continued to 
the complex w-plane so g{uj) can be written as 5(1^) — 
Tr-^[ — 7^ — ^T ,} , . J with poles u^ni 2 = idj^i'y). The 

2-772 La; — (u;— 27) lj — (Lj-t-27) -I ^ /^J-,^ \ // 

simplest case of model A, Eq. ([Ij, is when the external 
forcing is identical for each oscillator, h{I) = 5{I — e). 
This leads to trivial dynamics, solved by simply making 
the reference frame rotate at the TD frequency w + /(i), 
where Co is the mean of g(a;) and f + ft = ex. 

The non-autonomous (NA) dynamics arises for non- 
identical forcing. We first assume strengths proportional 
to frequencies, i.e. u}(t) — uj[l + exit)] with a constant 
e. The integration in Eq. Q is now only over w and, by 
closing the integral in any of the complex half-planes, is 
given by the residue of the encircled pole. As a require- 
ment from [18|], \a{uj,t)\ ^ as 3(w) — ?> Too, depending 
on which pole is encircled. The last limit transforms 
Eq. dH) into ^ = -Cj{t)a. Thus, for [1 + ex{t)] > 0, 
the encircling is around the pole Wp2 — {uj — 17), while 
for [1 + e2;(i)] < the upper half-plane encircling in- 
volves ujpi — (ilj + ij). The residue at these poles, 
z* = Q;(a) =F 17, t), is changed in Eq. ([7]) yielding 

f = -r[7|l + ex{t)\ + —{r^ - 1)], ^ = lo[1 + ea;(t)]. (10) 

The ansatz ([6]) holds only for nonidentical oscillators [19| , 
implying the requirement a)(i) ^ 0,Vi. 

Model A is also solvable with an independent 
Lorentzian distribution of forcing strengths. The fre- 
quencies follow Cj{t) = oj + Ix(t) and the mean and half- 
width of h{I) are / and 7/ respectively. The integrals in 
Eq. ([9]) can again be closed in the lower or upper complex 
half-plane, and the requirements for a(w, /, t) are similar 
to those in the previous case. Hence, the / integral for 
x{t) > is around the pole /pi = (/-!- 177) and around 
Ip2 — {I ^ ill) otherwise, while in the u integral the 
encircling is around the pole a;p2 = uj — i^. Thus, the 
residues give z* — a{ui — 17, / — 177, i), which is applied 
in Eq. ([7|), so we finally obtain 

r^ -r[7 + 77|a;(t)|-f — (r^-l)], i; ^ uj + ix{t).{\l) 

A similar analysis could be done for the other Lorentzian- 
like distributions of a; or / discussed in |18|. 
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FIG. 1. The time- varying mean field for model A, Eq. ([T]), 
resembles the externally applied cosine (a-c), or the chaotic 
forcing (d). Numerical simulations of the full system Eq. U]) 
(light grey) are in agreement with the low-dimensional dy- 
namics (dashed dark grey): Eq. (fTU)) - (a-b); Eq. (fTTj) - (c); 
and Eq. ^^ - (d). Adiabatic (dotted black), Eq. gH) - (a- 
b), and Eq. (|23[) - (c); and non-adiabatic evolutions (black), 
Eq. 1171) - (a-b), and Eq. ^^ - (c), confirm the limits of 
the reduced dynamics (see text for details). The distribu- 
tion of forcing strengths is: (a-b) same as the frequencies', 
K — 3.5, e = 0.6, fi = 5 and il = 0.03 respectively; (c) inde- 
pendent Lorentzian, K = 4.5, 7/ — 0.6 and n — 1; and (d) 
bimodal 6, K = 8, 'y = 1, -yi — 1 and 1 = 1. 



The only other analytically solvable form of model A 
that we are aware of is with multimodal (5-distributed 
external strengths. For simplicity we choose the bimodal 
function h{I) = ^[6(1 -I - 77) + 5{I - I + 77)]. The 
integral (jH]) now leads to 

^* = -[aiiuj-ijj -'-fi,t) + a2{Lb-ijJ + 11, t)], (12) 

with dynamics consistently described by the evolutions 
of ai^2 obtained from Eq. ([7]), 



dai^2 
dt 



= -{i[uj + (/ T li)x{t)] - 7}ai,2 + 

K 

+ — [ai + a2 - 01.2(0^1 + "2)*]- 



(13) 



For this case of model A, only, Choi et al. [22| made a 
bifurcation analysis near the limit rK <IC 1. 

Following the restrictions on x{t) in the problems an- 
alyzed in Fig. [11 we took x{t) = e cos fit for e < 1 in 
the case of strength proportional to the frequency, while 
for model A with independent Lorentzianly-distributed 
strengths, x{t) = 1-1- ecosfii and e < 1. Finally, for bi- 
modal (5-distributed strengths, the absence of restrictions 
on the external field allows it to be the x component of 
a Rossler oscillator [20J. In all the problems shown, the 
NA TD dynamics is revealed and fully described by the 
reduced NA low-dimensional system. A Runge-Kutta 4 
algorithm was used for numerical integration of Eqs. ([1] 
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FIG. 2. Time- varying mean field for model B, Eq. ((Sjl, fol- 
lows the external cosine (a) or the chaotic (b) forcing. Nu- 
merical simulation of Eq. ([2]) (light grey) coincides with the 
low-dimensional evolution (dashed dark grey), Eq. (|14p - (a) 
and Eq. HH) - (b). (a) Adiabatic (dotted black), Eq. ([2l|, and 
non-adiabatic evolution (black), Eq. (|20|) for constant forcing 
with K = 3, i^ — 0.1 and e — 0.33. (b) Bimodal (5-distributed 
strengths with K = 5, "fi — 1 and / = 0. 



El) over 100000 oscillators, with a time-step of 0.0025, 
while half- width and mean of the natural frequencies were 
■y — 1 and lu — 0, except where otherwise stated. 

We have also investigated the low-dimensional evolu- 
tion of NA model B, Eq. ([2]). Because all the couplings 
in the original model are equal, there is no qualitative 
difference between the situations with identical forcing 
to each coupling, and coupling-dependent forcing. We 
focus on the latter and proceed as for model A, yielding 

K 

f = -rij + —[1 + ex{t)]{r^ - 1)]; ^ = cb. (14) 

The analysis for multimodal (^-distributed strengths is 
very similar to that for model A ([T|). E.g. for bimodal 
h{I), the Eq. ([12]) holds again with ai,2 evolving as 

^ = -{ico - 7)ai,2 + \k[1 + (/ t li)x{t)] 

x[ai+ a2 — al2i(Xi+ 012)*]- (15) 

Nevertheless, for a Lorentzian distribution h(I), the re- 
quirement [18] that \a{uj,I,t)\ — )► as 3(/) -^ ±00 fails. 
In this limit Eq. ([8]) becomes: ^ = —^Icosm{za'^ — z*), 
and from there a{t) = ^e^"'^ ^ 0. This restricts us from 
applying ansatz ([5]) and the NA dynamics cannot be de- 
scribed. In contrast, the restrictions do not affect the 
NA parts of the other discussed variations of model B. 
To confirm this generality, x(t) for the problem shown 
in Fig. [2] (b) is a chaotic signal from a Rossler oscilla- 
tor. Similarly, the chosen amplitude of the cosine forcing 
in Fig. [2] (a) allows close-to-incoherent dynamics to be 
observed in some intervals, so that latter discussed limi- 
tations of the slow-fast approaches appear. 

A theorem in [19| states that Eqs. ^ [8|) asymptoti- 
cally capture all macroscopic behavior of the system as 
t -^ 00. Moreover the incoherent and partly synchronized 
states both belong to the manifold defined by Eqs. ([3 [8]) 
|18| . The initial incoherent state is set with uniformly 
distributed phases at time t — 0. Thereafter, the ansatz 
([71 [S]) and the evolutions (1101115^ describe our system con- 
tinuously, as confirmed by Figs. [T] and [21 



The most obvious feature of the results in Figs. [T] (a)- 
(c) and[2](a), is the low- frequency filtering of the external 
fields. E.g. the only difference between plots (a) and (b) 
of Fig. [l]is the frequency of the external forcing, while its 
influence is much more prominent in the latter. This is 
actually a well-known, but not much explored, character- 
istic of population models, and it is a direct consequence 
of their intrinsic transient dynamics [ij. In what follows, 
we adopt fast-slow reduction to simplify the evolution 
for simple periodic forcing. The reduction depends on 
the period of the external field T = 27r/i7, relative to the 
system's transition time, r, and has not been applied to 
similar systems. The exponential damping rate of the 
original system is defined by r [18| and t = l/\K/2 — 7|. 
For a system far from incoherence, K = 2j + 0(27), 
T ~ 1/0(7) holds, meaning that the transition time de- 
pends only on the width of the distribution of natural 
frequencies. Thereafter for this case, the system's re- 
sponse is adiabatic for slow external fields, $7^7, and 
non-adiabatic for fast, ^ ^ j. The dependence on 7 is 
removed by scaling the time and the couplings in the au- 
tonomous system, t = i/7, K = K/2j and t = 1/\K — 1\ 
(the scaled variables keep the same letters). 

For model A, Eq. ([TJ, with x{t) = e cos iU, after the ini- 
tial transition and in absence of bifurcations, the ampli- 
tude of the mean field consists of a constant term rg and a 
TD term Ar{t). For the non-adiabatic response, simula- 
tions, grey lines in Fig. [Ija), show that Ar{t) ~ l/H, and 
rg ^ Ar{t). Thereafter tq can be expressed as averaged 
over one period T = 27r/i7 of the oscillations of Ar{t). 
Thus, = -ro + K{4 -ro) [ll, or tq = ^/T^^TJk. 
Further, we apply r{t) w tq and ^ = -^ to Eq. ^ 
and then integrate it. From there Ar(i) = — 7'o^sinJ7i, 
and the magnitude of the NA response is 



A.ast = 21^1-1 (16) 

Hence the long-term non-adiabatic evolution follows 



rfast(<) = yi-;|(l-^sinffi). (17) 

The adiabatic behavior emerges through the introduction 
of a slow time-scale t' = fit, such that the system is 
constant on the fast time-scale t, and changes only in t' . 
Hence the l.h.s. of Eq. pO)) is zero, whence 



^slow \^) 



1 + e cos nt 



1- 



K 



(18) 



while, for the magnitude of the NA part, we obtain 



slow 



1 



1 -e 



1 



1 + e 



(19) 
K \ K ^ ' 

An analogous analysis can be performed for the appro- 
priate form of model B, Eq. ([2]), leading to the low- 
dimensional evolution for fast cosine forcing given by 



rfast(t) 



l + ^sinm)Jl-l, 



(20) 
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FIG. 3. (color online) Magnitude of the response, A, 

of the NA model A to the cosine forcing, Eq. ([!]). Exter- 
nal forcing strengths follow the distribution of frequencies, 
K = 4.25 and Q. G [10"^ 10^]. (a) Results from Eq. ^ 
for e e [0.05,0.99]. (b) Adiabatic (dotted black), Eq. ^ 
and non-adiabatic, Eq. (|19|l . (dashed black) evolution for 
e G {0.05,0.1055,0.2225,0.4693,0.99}, compared with the 
real dynamics (grey), Eq. (|10p . 



and for slow forcing 



^slow \^) 



•l- 



K{l + ecosnt)' 



(21) 



For the dynamics of model A with independent 
Lorentzian strengths and forcing, x{t) = 1 + I cos fit, 
the time is scaled by (7 + 7/), and the dynamics follows 



nast(0 



1 



7/ 



1 (1 

K^ f^(7 + 7/) 



sinr^i) (22) 



for fast forcing, while for slow driving 



^slow y^ ) 



I 1 7/ cos nt 

K K{j + 7/) 



(23) 



The adiabatic responses can also be obtained from the 
self-consistency of Eqs. (01 [S]) for stationary states of the 
mean field. Namely, assuming very slow dynamics of the 
external forcing, the system can be treated as quasista- 
tionary. This is similar to assuming stationarity on a 
fast time scale. Thus one obtains r — ^1 — 2j{t)/K{t), 
corresponding to the results (IT^ . (PT|) . 

Ah the evolutions for reduced dynamics, Fig. [Ha)-(c) 
and[2ja), are in line with the above analysis, confirming 
the interplay between external and internal time scales 
of the NA system. The magnitudes of the slow/fast re- 
sponses to cosine forcing are given in Fig. |3] for model 
A, Eq. ([T|), with forcing strengths following the distribu- 
tion frequencies. It confirms the obtained dependance of 
A on the frequency and amplitude of the external field. 
The low-frequency filtering mentioned before is also ob- 
vious. The transient behavior for slow and fast forcing 
can be seen in Fig. El^b), where A is shown for both the 
actual and the reduced dynamics. This plot perfectly 
matches the analytic limits for application of the reduc- 
tion approaches. Similar plots can also be obtained for 
the other problems analyzed. However, for coupling close 
to critical, the system's transition time increases and for 
K Ki Kc, T — >■ c». As a result, the slow dynamics fails, as 



shown in Fig. [U^a) when r is close to 0, unlike the case 
K = Kc + 0{Kc) given in Fig. [Ii;a)-(b) or Fig. [SJa) for 
r far from 0. 

With the analysis of the reduced dynamics, supple- 
menting the full low-dimensional description, all aspects 
of the TD KM have been demonstrated. The former is 
shown only for simple periodic forcing, but this does not 
decrease the generality of the reduction, since any exter- 
nal field can be represented by its Fourier components. 
Besides, these methods are of great importance in mod- 
eling systems with multiple time-scales of oscillation and 
interaction, such as e.g. the human cardiovascular system 
2J], or inhibitory neurons in the cortex j25| . 

In summary, we have characterised a new dynamics 
of interacting oscillators subject to continuous, deter- 
ministic perturbation. It consists of the dynamics of an 
external system superimposed on the original collective 
rhythm and was missing from earlier models [3| , possibly 
leading to wrong interpretation of some real dynamical 
systems. We have derived the impact of the forcing and 
evaluated the effect of its dynamics, amplitude and dis- 
tribution. Thus, we have proposed a generalization of the 
Kuramoto model that encompasses NA systems [2a] and 
is directly applicable to any thermodynamically open sys- 
tem. E.g. the observed time- variations of brain dynamics 
can be easily explained as a consequence of TD frequen- 
cies or couplings of the single neurons, where the source 
of the external variation could be due to anaesthesia [5|| , 
event- related [7|, or due to some influence from another 
part of the brain. In particular, the stable, time- varying 
mean field can now be reconstructed and, in this way, a 
large range of systems tackled by the Kuramoto model - 
spanning from a single cell up to the level of brain dy- 
namics - can be described more realistically. 
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